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Abstract 

We present an algorithm which produces a decomposition of a regular cellular com- 
plex with a discrete Morse function analogous to the Morse-Smale decomposition of 
a smooth manifold with respect to a smooth Morse function. The advantage of our 
algorithm compared to similar existing results is that it works, at least theoretically, 
in any dimension. Practically, there are dimensional restrictions due to the size of 
cellular complexes of higher dimensions, though. We prove that the algorithm is 
correct in the sense that it always produces a decomposition into descending and 
ascending regions of the critical cells in a finite number of steps, and that, after a 
finite number of subdivisions, all the regions are topological discs. The efficiency of 
the algorithm is discussed and its performance on several examples is demonstrated. 
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1 Introduction and motivation 



Given a smooth manifold M witli a Morse function F defined on it, classical 
Morse theory [1], [2] is a powerful tool for investigating the behavior of the 
function F, as well as the topological properties of M. The critical points of 
F together with their indices determine a handlebody decomposition of the 
domain M. A different decomposition of M, carrying topological as well as 
geometric information, is obtained from the stable and unstable manifolds 
of the critical points. The intersections of these are regions where the func- 
tion behavior is uniform in the sense that all gradient paths have the same 
asymptotic behavior. If the function is Morse-Smale, that is, if the stable and 
unstable manifolds of the critical points intersect transversely, then these re- 
gions form the Morse-Smale CW-complex. Such a decomposition enables a 
thorough understanding of the domain M, as well as of the gradient flow of 
the function F. 

Consider for example the function shown on figure [H Assume that the function 
models a geographic terrain, and that our task is to find a path starting at 
a point (x, y) in the descending disk of one of the maxima and ending in our 
preferred maximum (xo,yo)) which is traced optimally with respect to some 
criterion, for example height variation. This can be reconstructed from the 
descending and ascending disks, but these are computationally expensive. The 
problem is even more difficult if, instead of the function /(x, y), only a sample 
of points on the surface is given. What we need, is a discrete approximation of 
the ascending and descending disks of F. Once this is given, the required path 
can be constructed so that it first reaches a saddle in the common boundary 
of both maxima, and then ascends towards the preferred one. 

In this paper we propose an algorithm based on the discrete Morse theory of 




Figure 1. A model of a geographical terrain 
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Forman [3], [4] to solve this problem. Discrete Morse theory is a PL analogue 
of classical smooth Morse theory which has gained wide popularity and is 
used in topological data analysis [5] , [6] , [7] , [8] as well as in addressing purely 
theoretical topological and combinatorial problems, as for example in [9], [10], 
[11]. 

In order to use discrete Morse theory for analyzing data, the initial data are 
first extended to a discrete Morse function on a triangulation (or, more gen- 
erally, regular cellular decomposition) of the domain. This can be done using 
for example the algorithm of [12] (which is used in our implementations). The 
algorithm of this paper provides a further step by constructing the descend- 
ing and ascending regions of the discrete Morse function and, from these, the 
discrete Morse-Smale complex. 

The construction of discrete descending regions of a critical cell is motivated 
by the definition of descending disks in the smooth case: as unions of V- 
paths, which are the discrete analogue of gradient paths, starting in the critical 
cell. Discrete ascending regions are constructed using the same procedure on 
the dual l/-paths in the dual complex K*. The resulting decomposition into 
descending regions is similar to the smooth case in top dimension where the 
obtained regions are disjoint topological disks. In lower dimensions they are 
not necessarily disjoint, and not always disks, since, unlike gradient paths 
in the smooth case, V-paths from different critical cells or even from the 
same critical cell ending in the boundary of a specific critical cell s can merge 
before reaching this boundary. We show, though, that after a finite number of 
subdivisions such merges can be eliminated, so that, as in the smooth case, 
the descending regions and the ascending regions form two families of disjoint 
topological disks. The regions obtained typically do not intersect transversely. 
This is not a major problem in the discrete case, though. The intersections 
nevertheless provide a discrete Morse-Smale decomposition of the cells of the 
domain into regions where, like in the continuous case, all y-paths come from 
and tend towards the same two critical cells. 

In [11], Theorem 3.1 the connection between discrete and smooth Morse the- 
ory is described in detail. It is shown that a smooth Morse function F on 
a closed Riemannian manifold M can be approximated by a discrete Morse 
function / on a triangulation of M so that critical points of F correspond 
to critical cells of /. If, in addition, the function F is Morse-Smale, the V-paths 
connecting any pair of critical cells p and q such that dim(p) = dim(g) + 1 are 
in bijection with integral curves of the gradient vector field of F between the 
corresponding critical points. 

Several other algorithms for constructing discrete ascending and descending 
disks, as well as a discrete Morse-Smale complex exist. A classical example 
of an application of such a decomposition is the watershed segmentation al- 
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gorithm in digital image analysis [13,14]. In [15] a discrete Morse-Smale de- 
composition of a surface was constructed and applied to geographic surface 
modeling, in [16], using a different approach, such a decomposition was used in 
molecular modelling. In [17] 3D Morse-Smale complexes are used in volumet- 
ric data analysis. In [18] an algorithm for constructing a quasi Morse-Smale 
complex, a combinatorial analogue of a 3D Morse-Smale complex, is given and 
in [19] an extension of this algorithm is applied to several well known data 
sets. An overview of computational Morse theory with a large number of addi- 
tional references and applications can be found in [20]. All existing algorithms 
are restricted to 2 or 3 independent variables, though, while real-world data 
often depends on more than 3 variables. 

Our algorithm works, theoretically, in any dimension, and we have implemen- 
tations in OCAML, C-I--I- and C# that accept simplicial or, more generally, 
polyhedral complexes without dimensional restrictions as input. In practice, 
however, the implementations can be used for simplicial complexes of dimen- 
sion 6 or less, since the number of cells in a complex grows exponentially with 
the dimension, and the complexity of the algorithm depends on the number of 
cells. One of the applications described in the last section uses a 4-dimensional 
data set. 

An additional advantage of discrete Morse theory, in comparison to other 
approaches, is an efficient and elegant mechanism for dealing with noise which 
produces, after canceling pairs of critical cells, a simplified function in the 
sense that the number of critical elements is reduced. This is similar to handle 
sliding in computational Morse theory of [18] and [15], and allows the use of 
persistence [21], [22] which has proven to be extremely useful in topological 
analysis of data [23], [24], [19]. 

The paper is organized as follows. In Section [2] we give a short overview of 
discrete Morse theory and introduce notation. In Section [3] we describe the 
construction of descending and ascending regions of critical cells, and discuss 
the complexity and the properties of the resulting decomposition. Finally, Sec- 
tion m contains several demonstrations and applications of our algorithm. We 
first give a short description of an application of our algorithm to qualitative 
data analysis in artificial intelligence presented in [25]. Next, an analysis of 
the response of a mechanical system consisting of a cart and a rod attached 
to its top is described. Our third example concerns the problem, presented in 
the beginning of this section. A procedure for constructing an optimal path 
(with respect to some optimization criterion) using the decomposition into 
descending regions is described, which is well defined also in the case, where 
the descending and ascending regions are not disks. Finally, we use our de- 
composition to construct a macroeconomic model on a 4-dimensional data 
set. This is part of an ongoing project with the goal to find and algorithm 
for controlling the key parameters of the model which contribute to long term 



4 



economic growth of a country. 



Part of this work was the resuh of a cooperation with the Laboratory for 
artificial intelhgence at the Faculty of Computer and Information Science, 
University of Ljubljana. We would like to thank the members of this lab for 
the permission to use data gathered in their experiments. Our thanks also 
go to the anonymous referees for valuable comments which have helped to 
improve this paper. 



2 Discrete Morse theory 

A discrete Morse function on a regular cellular complex K associates to every 
cell a & K a. real number /(cr) such that in most cases / increases with 
increasing dimension except, possibly, in one direction. More precisely, for 
every cell r^^^ G K of dimension p the number of its co dimension 1 faces 
with values of / greater than /(cr) is at most one, and also the number of its 
codimension 1 cofaces with values of / smaller than /(ex) is at most 1. That 
is: 

6(r) = #{z/(^-i)|z/<r,/(z/)>/(r)}<l 
a(r) = #{^(^'+1) I r < a, /(r) > f{a)} < 1. 

The numbers a(r) and 6(r) can not both be one. Since, if there exists a face 
v < T such that /(z/) > /(r) as well as a coface a > t such that /(a) < /(r), 
then for any t' < a such that z/ < r', as in figure [2] (where the arrows indicate 
the direction of function descent) it follows that 

f{r') > f{v) > fir) > f{a) > f{r') 

which is a contradiction. 




Figure 2. A discrete Morse function can not have both a(r) and 6(r) equal to 1 

Because of this, the cells of K split into three subsets K = AVJ B U C {m 
the terminology of [12]). If both a(r) = 6(r) = then r is a critical cell of 
index p = dimr, and C is the set of all critical cells. If not, either a(r) = 1 
or 6(r) = 1, and r is a regular cell. A regular cell r belongs to the set A if 
a(r) = 1 and to B if 6(r) = 1. 
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For each t & A there exists precisely one a E B such that r < a and /(r) > 
/((j). The map V : A ^ B which associates t^^^ to a^^^^^ is called a discrete 
gradient vector field, and points in the direction of steepest descent from the 
cell r. 

It is convenient to imagine a discrete gradient vector field not as a map but 
as a collection of pairs (T^^\a^^^^^) such that V{t) = a, or arrows from r 
to a, in the direction of function descent. This representation of a discrete 
gradient vector field is used in all figures in this paper. The cells of K which 
are unpaired are precisely the critical cells. 

In classical Morse theory the gradient vector field points in the direction of 
steepest ascent of the function /, and its gradient paths (traced in the direction 
opposite to the gradient flow) are the paths of steepest descent. The discrete 
analogue is a y-path, which is a sequence of cells 

such that V{Ti) — Ui and Tj+i is a face of CTj different from for i e {1, .., r}. 

As in the classical case, function values decrease along a F-path. This implies 
that a l^-path of a discrete gradient vector field V arising from a discrete 
Morse function can not be closed, that is, the discrete gradient vector field can 
not contain any cycles. Forman showed that this is characteristic for discrete 
gradient vector fields: a pairing V arises from a discrete gradient vector field 
of a discrete Morse function on K if and only if it contains no nontrivial closed 
y-paths. 

In classical Morse theory, two critical points pi and p2 of a Morse function F 
on M of indices differing by 1 can be canceled when there exists exactly one 
gradient path between them. That is, using handle sliding, F can be smoothly 
deformed to a function F' which has exactly the same critical points as F 
except pi and p2 that become regular. This is true also in the discrete case. If 
the critical cells a and r are connected by exactly one V-path starting in any 
cell of da and ending in r, i.e. 

then the pair cr, r can be cancelled. This is done simply by modifying the 
discrete gradient vector field V along the connecting paths to 

for i G 0,...,r + l. This procedure can also be described as switching the 
direction of the arrows along the connecting V-path. Since there are no other 
V-paths connecting Tq^^ < da and r*^P^ this can not create a non-trivial closed 
y'-path and so, according to Forman's characterization, the modified V is the 
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discrete gradient vector field of some Morse function for which the cells r and 
a are no longer critical. 



3 Algorithm 

In this section we present our algorithm which, on the basis of a discrete Morse 
function /, produces two decompositions of a regular cellular complex K with 
\K\ a manifold. First, a decomposition of K into discrete descending regions 
of the critical cells of / is constructed. Each of the regions contains exactly 
one critical cell s*-^-*, the maximum of the region, and a collection of regular 
cells of dimension less than or equal to p. 

In the second step we decompose, using the same procedure, the dual complex 
K* into discrete descending regions of the function — /. This gives a decom- 
position of K into discrete ascending regions of the critical cells. 

3.1 Constructing the descending regions 

In the smooth case, the unstable manifold of a critical point c of F is the union 
of all gradient paths of the flow — grad F starting close to the point c. In the 
discrete case the direction of function descent is indicated by the discrete 
gradient vector field V , which plays a similar role as the gradient vector field 
in the smooth case. The discrete descending region D{s^^'>) of a critical cell 
s^P'^ thus contains all regular cells which appear in a V^-path starting in the 
boundary of s^'p\ The descending region is constructed in two steps. In the first 
step the frame which contains all cells of maximal dimension is constructed, 
and in the second step the cells of lower dimension are added. 

The critical cells are first ordered by ascending dimension so that, when pro- 
cessing a critical cell of dimension the descending regions of all critical cells 
of dimension less than d have already been constructed. This will be very 
important in the second step of the algorithm. 

Let s'-^-' be a critical cell. The frame f r(s) consists of all regular p and {p — 1) 
dimensional cells appearing in a V^-path starting in the boundary of s. This 
part of the algorithm is a simple breadth-first search through the regular 
{p — l)-faces of the included p-cells which terminates when no more regular 
(p— l)-cells belonging to the set A (that is, no more {p— l)-dimensional arrow 
tails) can be found. Since a cell in the frame can typically be reached by several 
paths in V , a set structure is used for storing the cells to eliminate duplicates. 

This step of the algorithm is illustrated in figure O On the left, the discrete 
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gradient vector field is shown, the descending frame of the maximum is in the 
middle and the descending frame of the saddle is on the right. 
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Figure 3. Building the frame of descending regions in a cellular complex 

The frame f r(s'''''-') contains only p and p — 1 dimensional cells. The next step 
is to decide which additional faces of the cells in f r(s) should be included. 

Consider a y-pair {a, (3) incident to some simplex from fr(s^*''*). Since V{a) = 
P, function values decrease from a towards p. This implies that P belongs to 
the descending region D{s) when a does, so the pair (a,/?) is included into or 
excluded from the descending region together. It is included if all codimension 
one cofaces of a except P have already been included. This is implemented 
by the following recursively used routine. Let S be the set of all {p + 1)- 
dimensional cofaces of a, different from P, and incident to some p-dimensional 
cell from fr(s). The routine checks if all elements of S belong to D{s). It ends 
when either a cell 7^'^+^) G S" is found that belongs to the descending region 
of some other critical cell of dimension less than p, in this case the pair 
is not included, or all cells of 5* have been included in D{s), and in this case 
the pair {a,P) is also included. 

On figure H] the left picture again shows the discrete gradient vector field. In 
the middle picture, part of the descending region of the maximum is already 
constructed (the 2-cells forming the frame are shaded, and the 1-cells that 
are already included are bold), and the pair (a,/?) that is processed is the 
diamond-shaped point and the dashed edge. Since all edges that have this 
point as a face are already included in the descending disk, the pair {a,P) 
is added. On the right, the situation is different. The descending disk of the 
saddle has already been constructed and contains the dashed edge p. Since P 
already belongs to the descending region of the saddle, the pair (a, P) will not 
be included. 

Proposition 1 The construction of the descending region of a critical cell 
ends in a finite number of steps. 



PROOF. During frame construction we simply follow the direction of the 
discrete gradient vector field V. Since V does not contain nontrivial closed 
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Figure 4. Completing the descending region 

paths and the complex K is finite the frame is completed in a finite number 
of steps. 

The second step of the construction, where lower dimensional cells are added, 
is recursive. The algorithm runs through all y-pairs in the boundary of frame 
cells. For each pair {a^'^\ P^'^'^^^) which has not been processed it checks if all 
cofaces of a of dimension {d+1) are already included in the descending region. 
If it encounters a coface which has not been processed it checks that first. This 
implies that the algorithm always proceeds in the direction opposite to V, so 
no cycles can be created. 

Since the number of pairs is finite, the algorithm ends in a finite number 
of steps. Notice that the algorithm can start processing pairs in any order, 
independent of dimension. 

3.2 The boundary 

The algorithm described in the previous section constructs discrete descending 
regions of a discrete Morse function on a complex K with \K\ a manifold with- 
out boundary. For manifolds with boundary, the boundary must be processed 
separately. 

Denote the discrete gradient vector field of the restriction of / to the boundary 
by Vqk- a critical cell u of Vqk is 

(1) either a critical cell in V 

(2) or is paired to a higher dimensional cell that does not belong to the dK; 
such cells will be called boundary critical cells. 

In the first case no extra processing is needed, as is also a critical cell in V 
and its descending region has been constructed in the previous step. 

In the second case, the discrete descending region of a boundary critical cell 
ly^P^ is constructed in two stages. First, the discrete descending region of the 
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critical cell dK is found using the algorithm of the previous section. In 

the second stage, all cells of dimension p in the discrete descending region of 
dK^ that are paired by to a higher dimensional cell (3 in the interior, 
are considered. The cell {3 is processed as a critical cell and its descending 
region is constructed. The union of all regions obtained in this way constitutes 
the descending region of the boundary critical cell v. 

For example, consider /(x, y) = a; + y on a triangulation of the square [0, 1] x 
[0, 1]. The discrete gradient vector field contains all cells except a minimum 
on the boundary which is the only critical cell 




O 



Figure 5. Discrete gradient vector field on the square and its boundary 

The discrete gradient vector fields of / and / 1 qk are shown on Figure [51 The 
maximum /x of /| ax is a boundary critical 1-cell. Its descending region in OK 
is all of OK except for the minimum. A boundary cell a which is paired with 
an interior cell (3 in V is included into the descending region of /i together 
with (3 and its descending region. Finally, the descending region of /i is the 
entire square without the minimum. 

3.3 Ascending regions 

The ascending regions of the critical cells of / on are obtained from the 
dual complex K* with the dual the discrete gradient vector field V*. 

Viewing K as an abstract regular cell complex determined by its cells with 
specified dimension and the face relationship between them, the dual complex 
K* is obtained simply by reversing the dimensions and the face relationship 
on the cells. That is, a cell a*^^^ G K corresponds to a cell a*^*^"^^ G K* ^ where 
n is the dimension of K, and if a is a face of (3 in K, then a* contains [3* 
as its face in K*. In the geometric realization \K\ of K, the vertices of K* 
correspond to barycenters of the n-cells of K, while a cell of higher dimension 
is spanned by a collection of vertices of K* such that the corresponding n-cells 
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of K have a nonempty intersection. If l-ft'l is a PL-manifold without boundary, 
then \K*\ = \K\. If \K\ is a PL-manifold with boundary, then K* is not a 
cell decomposition of any more, since the cells of K* which are dual to 
boundary cells of K do not have a sufficient number of faces (for example, a 
1-cell of this type has only 1 face). For our purpose K* is just a combinatorial 
object, though, so this geometric property does not represent a problem. 

Also the dual discrete vector field V* is obtained simply by reversing the 
arrows in V. That is, two cells a* and (3* are paired in V* if the cells (3 and a 
are paired in the original discrete gradient vector field V. Thus, a cell a in K 
is critical of index p if and only of the cell a* is critical of index {n —p) in K*. 

The discrete ascending region of a critical cell a of is the dual of the discrete 
descending region of the critical cell a* of K* (see figure [6]). 
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Figure 6. Ascending region of a saddle 

3.4 Consistency 

The descending and ascending disks of a smooth Morse function / on a man- 
ifold M form a decomposition of M into topological disks. If the function is 
Morse-Smale, their intersections form the Morse-Smale complex on M. In the 
discrete case, the descending and ascending regions obtained are always disks 
only in the top dimension, while in lower dimensions they might not be disks. 
In addition, they might not be disjoint, that is, a cell can belong to the de- 
scending or ascending disk of more than one critical cell. We first prove that 
the descending regions (and therefore also the ascending regions) of all critical 
cells form a covering of the underlying PL manifold. 

Proposition 2 Let K be a regular n- dimensional cellular complex such that 
\K\ is a manifold without boundary and let V be the discrete gradient vector 
field of a discrete Morse function on K . Then every regular cell a & K is 
contained in the discrete descending region of a critical cell (3. 

This is true also for manifolds with boundary, if the discrete gradient vector 
field of f on K is an extension of a discrete gradient vector field of f \qk- 
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PROOF. We first prove that every regular top dimensional cell a'^"^ is in- 
cluded in the frame of some top dimensional critical cell. A regular cell of 
top dimension is paired either to a cell 70" which is a face of at least one 
n-dimensional cell cci other than ckq, or it is paired to a boundary critical 
cell of top dimension in dK. In the second case, q;^") belongs to the 

descending disk of z/*^"^^\ In the first case, if ai is critical, then ao lies in 
its descending region, and if it is regular, we repeat the same process for ai. 
Since the discrete gradient vector field V does not have cycles and K is a. finite 
cellular complex this process must end in a finite number of steps at a critical 
cell CKj-"^ or a boundary critical cell, and ao belongs to its descending disk. 

Assume now that descending disks of all critical cells have been built and that 
there exists a regular pair {a^\ that is not included in any of them. 

Since every cell is a face of some top- dimensional cell, a was processed dur- 
ing the construction of the descending region of some cell a. Since it was not 
included, there exists a coface 71^^^^ of a, such that 71 is included in the de- 
scending region D[ai^^^), where rii < n. In the construction of the descending 
region D{a^^'^) also a was processed since it is also a face of 70, and not in- 
cluded. So there must exist a coface 72^'''^'' of a which belongs to some other 
descending region D((T2"^^), where n2 < ni. This argument can be repeated 
indefinitely, and we obtain an infinite sequence of critical cells ctq, cJi, a2 . . . 
with strictly decreasing positive dimensions which is a contradiction. So every 
regular pair {a'^\P^~^^^) is included in the descending region of some critical 
cell. 



In smooth Morse theory the stable and unstable manifolds of a critical point 
are constructed from gradient paths, that is, integral curves of the gradient 
vector field, which begin or end at this point, respectively. An important prop- 
erty of integral curves of a smooth vector field is that they are pairwise disjoint 
- two different integral curves can be arbitrarily close but do not merge. The 
discrete analogue of a gradient path is a l^-path, and since there is no such 
thing as 'arbitrarily close' in the discrete world, a \/-path typically splits into 
multiple paths, and several V-paths can merge into one path. Due to this, dis- 
crete descending and ascending regions of different critical cells might not be 
disjoint, since y-paths beginning in different critical cells might merge. They 
also might not be topological disks, except in the highest dimension. 

Proposition 3 The descending region of a critical cell of maximal dimension 
(i.e. a maximum) collapses to the maximal cell. The same is true for ascending 
regions of a critical cell of dimension (i.e. a minimum). 
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PROOF. For every cell a of maximal dimension, its pair is a cell of codi- 
mension 1 which is either a boundary critical cell or it is the face of precisely 
one other cell of maximal dimension. Because of this there exists precisely one 
V^-path starting at some (possibly boundary) critical cell s of maximal index 
(i.e. a maximum) leading to a. Since a discrete vector field has no cycles, ev- 
ery such ^-path leading to a cell in the boundary of the descending region 
of s determines a sequence of elementary collapses starting in the boundary 
of the descending region and ending in s. A finite number of such \^-paths 
determines the necessary elementary collapses required to collapse the whole 
descending region to the cell s. 

For descending regions of critical cells of index smaller than maximal, i.e. 
saddles, this is not necessarily true, since F-paths not containing cells of max- 
imal dimension can merge. Figure [7] shows a section of a 3-dimensional cubical 
complex where V^-paths in the descending region of a saddle split and merge 
again, so that the union is not a disk. 




Figure 7. Discrete descending region of a 2-saddle with splitting and merging 
y-paths 

It is possible, however, to modify the decomposition of the cellular complex 
and the discrete gradient vector field so that the merging point of descending 
paths is pushed further towards the boundary of the descending region, while 
the discrete vector field outside the open star of the merging point is left 
unaffected. The idea is similar to sphtting multiple saddles in [15]. 

The left example on figure [8] shows a section of a 2-dimensional simplicial 
complex where the descending region of a critical cell of index 1 is topologically 
a 1-sphere. The right side of this figure shows the modification which pushes 
the merging point further along descending values of /. The figure shows only 
one step of the pushing procedure, which can be applied as many times as 
necessary. 

Here is a description of the required modification of a regular cellular complex 
in the general case. Assume that two l^-paths 
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. . . 7j(^-^) < > rf^-'^ = r < a^^^ > t^J^'^ • • • 

merge in the cell t^p~^^ which is paired with a^P\ Let 5* be star(r) and let pi 
and p2 be two disjoint paths in S from the cell cTi_i to the cells a and aj_i 
respectively. 

The required modification consists of replacing t^p~^^ by two parallel copies 
T^P^^^ and t'^p~^\ and cr*-^^ by two copies a^^^ and cr'*-^-' with da' containing the 
same cells as da, except for r which is replaced by r'. In addition, two new 
cells are added: a cell u^^^ which has cells t^p~^^ and its boundary, 

and a cell h^p~^^^ with u, a and a' in its boundary. We modify the star S in the 
following way 

• in the boundary of all p-cells in the path pi (except a) the cell r is replaced 
by r', 

• in the boundary of the last {p + l)-cell in the path pi a is replaced by a', 

• z/ is added to the boundary of the first (p + l)-cell in the path p2, 

• r' and u are added to the boundary of all other cells in 5*. 

The modified field V coincides with V on both copies of r, a and their bound- 
aries, and it pairs u with fi. The new field V now has two \^'-paths, which 
possibly merge in some cell Xj+fc = rj_,_; (or are disjoint). So we have pushed 
the merging point further along the discrete gradient vector field V. Also note 
that this step does not affect ^-paths originating from critical cells in dimen- 
sion p — 1 or lower and does not create any additional merge points for critical 
paths originating from critical cells of dimension p. 

Proposition 4 Let K be a finite cellular complex without boundary and let 
V be a discrete gradient vector field on K . Then there exists a finite sequence 
of steps described above which modify K into K' and V into V in such a 
way that the set of critical cells remains unchanged, and that the descending 
regions of all critical cells with respect to V are topological disks. 



PROOF. The proof goes by ascending dimension of the cells of the complex. 
Since there is nothing to do for 0-dimensional critical cells, we start with 
critical cells of dimension 1. Let rf^^ be the cell where \^-paths originating 
from a critical 1-cell s merge and let D be the descending region of s. Since K 
is finite, a finite sequence of steps described above pushes the cell r^°^ into the 
boundary of D. We repeat this process for all such merges in the descending 
region of s. Let K" be the resulting complex, V" the resulting gradient vector 
field and D" the discrete descending region of s in K' . Since D" does not 
contain its boundary, and since all merging points have been pushed into the 
boundary we have pushed the cells where V^"-paths originating from s merge 
away from D" . Now the same argument as in the proposition [3] can be used 
to show that D" collapses to the critical cell s and is therefore a topological 
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disk. Also, all K-paths coming from critical cells different from s have been 
pushed off D", so D" does not intersect any other descending regions. 

Since all modifications take place in the interior of the descending disk of 
s, all other descending regions remain unchanged. And since no additional 

merge points have been created in the process it follows that, after repeating 
this procedure for all critical 1-cells, all descending regions originating from 
critical cells of dimension one are topological disks. 

Precisely the same arguments work for critical cells of higher dimension. As- 
suming that descending regions for critical cells of dimension less than p have 
already been processed and are therefore topological disks, the descending 
regions for each critical cells of dimension p is processed in the same way as 
above. Since all merging points are pushed to the boundary and since no other 
descending regions are affected, all descending regions of critical cells of di- 
mension p are transformed one by one into disks. Moreover, since the pushing 
step does not affect the descending regions of critical cells of dimension less 
than p, descending regions of all critical cells of dimension up to p are now 
topological disks. 





1 * O 












V 


T 






t' 





Figure 8. One step of pushing a merging point towards the boundary 



3.5 Complexity 



Let m be the total number of cells of fC, n the dimension of the complex i^, 
the number of cells of dimension d, the average number of codimension 
1 faces (if X is a simplicial complex = d -|- 1), p,^ the average number 
of codimension 1 cofaces and Pmax = niax(po,Pi, ■ ■ ■ Let Cj be the total 
number of critical cells of the discrete Morse function / on and q the 
number of critical cells of dimension d. 

The frame of a d-dimensional critical cell s^'^^ consists of d and [d — 1)- 
dimensional cells that belong to F-paths starting in the boundary of s. For 
each (i-dimensional cell in the frame, the set of its (d — 1) -dimensional faces 
has to be found, and for each face its pair in V has to be found. Since faces 
and their pairs can be found in constant time, this takes at most a x r^^ x 
operations for a suitable constant a. 
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Building all frames of critical cells thus takes Ya=i a x rrii x ri operations. In 
the case of simplicial complexes, the number of operations is at most (n+ 1) x 
a X m X n. The complexity of this step is therefore x m). 

To complete the descending regions requires determining whether a regular 
cell q that is incident to the region of some descending disk is included in 
this region. First, its V-pair p has to be found which can be done in constant 
time. Next, a list of all co-faces of the lower dimensional cell of the pair is 
required which can also be found in constant time. Each cell in this list is 
then processed in the same manner. On average it takes pj operations for each 
cell of dimension j, so we need at most Ylf^i a" xpiXuii operations to complete 
the descending region of one critical cell of dimension d, where a" is a constant. 
Since we can not use previous results when completing frames for other critical 
cells some cells are typically checked more than once. So completing a frame 
for all critical cells takes at most Cf Z)iLi o;" x Pmax x m — Cf x Pmax x m x n 
operations. 

For simplicial complexes the complexity of the algorithm is altogether 0{nm x 
(n + CtPmax))- Note that the algorithm complexity increases exponentially with 
the dimension of the complex, since the number of cells grows exponentially. In 
the special case where K is the Delaunay triangulation on a set of points in M", 
an upper bound for the total number of cells is 0(^1^"/^^) where v is the number 
of points. In this case the algorithm has complexity 0{v^'^^'^^n x (n + CtPmax))- 



4 Examples 

The algorithm presented in this paper has been applied to data sets from 
different domains. In this section we present some of them. 

4.1 The QING algorithm 

We first mention the algorithm QING which is an application of our algorithm 
to AI from [25]. The algorithm uses discrete Morse theory to reconstruct the 
critical cells of a discrete Morse function obtained from the values of a sampled 
function. The descending disks obtained using an early implementation of 
this algorithm were used for the construction of a qualitative graph, that is, 
an undirected graph G — {Vq, Eq) that represents connections between the 
critical cells. The set of vertices Vg is the set of all critical cells. Two points a 
and j3 from Vg with dim a > dim f3 are connected if, f3 is in the boundary of the 
descending region of a. This graph was then used in learning in a qualitative 
model. 
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The QING algorithm, together with the parametric discrete Morse theory 
of [6] was used in a learning scheme presented in [26] designed to teach a 
robot the concept of occlusion. 

4-2 Mechanical system response 

Our next example is obtained from a model constructed for the purpose of 
qualitative modeling of data in artificial intelligence. The model is taken from 
[27]. We would like to thank the authors for permission to use the data. The 
data represents measurements obtained from a system consisting of a cart and 
a rod attached to its top in such a way that it can fall either backwards or 
forwards. The cart is set on ice (or other low-friction surface) and the rod is 
positioned upwards. As the rod begins to fall, the cart responds, accelerating 
in the opposite direction of the fall. 

Our discrete Morse function models the response of the cart in terms of the 
acceleration x, depending on the rod angle if and the angular velocity (p. The 
measured values of x are given on a square grid. A discrete gradient vector field 
was first constructed on a triangulation of the grid using the algorithm of [12] 
producing, after cancellation, three maxima, four saddles and two minima. 

On figure [9] the discrete descending regions of all three maxima are on the 
left, and discrete descending regions of the four saddles are on the right. The 
descending regions form a decomposition of the triangulated area into disjoint 
discs. Note that the critical cell on the right side of the image is a boundary 
saddle. Its descending area is correctly generated and separates the descending 
regions of the two maxima. 

The exact equation describing this model is 



Figure fTOl shows the graph of this function (with suitable values of a, 6, c and 
d) with the descending discs mapped onto it. 

The motivation for this example comes from qualitative modeling and simula- 
tion in artificial intelligence [28], as a step further after the algorithm QING. 
The idea is, that understanding the response of the cart can provide an opti- 
mal strategy for controlling the system in the unstable equilibrium. Given a 
state of the system described by values of ip and 0, a V^-path from the given 
state towards this point encodes a strategy for control. Depending on some 
optimization criterion, an optimal strategy can be chosen. In the next example 
we give a brief explanation of a controlling strategy towards a maximum. 



aLfP' sm{ip) — 6sin(2(y9) 




X = 



c — d cos^ (f 
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Figure 9. Discrete descending regions of maxima and saddles of the cart function 



Figure 10. Graph of the cart function with a decomposition into descending disks 
4-3 Optimal path construction 

In this example, a discrete approximation of the Morse-Smale decomposition 
of the smooth function on figure [T]is first constructed. The function is sampled 
on a random set consisting of 10^ points in the domain. The Delaunay trian- 
gulation on these points consists of approximately 4* 10^ simplices. A discrete 
Morse function on this triangulation is constructed. The discrete Morse-Smale 
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complex is shown on figure [TTJ The required time to build the Morse-Smale 
complex on a PC was a few seconds. The construction of the discrete Morse 
function was computationally the most expensive and required (on a laptop) 
approximately 10 minutes. 

We also present a procedure for designing an optimal path, with respect to 
an optimization criterion involving height variation, from a given a point A = 
[x, y) in the domain, ending at a preferred maximal cell (3 (for example at 
the highest one). A slight modification of the described procedure can be 
used for constructing a path towards any critical point or even towards any 
point in the domain. Note that the procedure described can be applied in any 
dimension and can be therefore used to solve similar control problems in a 
more general setting. Also note, that the path constructed does not have the 
smallest possible height variation. Such a path would involve the construction 
of level curves which is not the object of this paper. 

Assuming that a triangulation of the domain, a discrete Morse function on it, 
and a discrete Morse-Smale decomposition are given, we first determine the 
set of critical cells S with the property, that A lies in their descending regions. 
Since every cell is contained in at least one descending region, the set 5* is 
nonempty, and because A may belong to several descending regions, the set 
S may contain more than one element. 

Let a be a cell of the same dimension as (3 with a & a. If (3 & S then a is a 
cell of maximal dimension in the descending disk of (3, so there exists a unique 
l^-path from (3 to a, and we simply climb up from a along this path. 

If (3 ^ S, we construct the qualitative graph, as in paragraph 14. 1[ The next 
step is to search through paths in the graph G starting from any point 7 G 5* 
and ending in the required destination (3 to find the path s which minimizes 
height variation. Note that any other criterion can be used to choose the 
optimal path at this step. 

The required optimal path consists of two parts. The first part connects a 
with the closest saddle 7 of index less than /? in s which lies in the boundary 
of the descending region of the initial point of s. The second part of the path 
follows s from 7 to /3. 

The cell a is connected to 7 in the following way. First, all \^-paths from a are 
constructed. If any one of them intersects the ascending or descending disk 
of 7, we follow this path to the boundary, and then ascend or descend to 7. 
If none of them do, we make a step along the unique V^-path from the initial 
point of s to a, and repeat the procedure. After finitely many steps we either 
reach the descending or ascending region of the initial point of s, or a, from 
where a V^-path to this point exists. 
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Figure 11. 3D image of the terrain function on figured] 
4-4 A macroeconomic model 

In our last application, macroeconomic financial indicators of European coun- 
tries are analyzed. The data was retrieved from the publicly available database 
Eurostat [29]. This example is part of an ongoing project involving a quali- 
tative analysis of the effect of key macroeconomic indicators on long term 
economic growth of a country. The motivation for this project is the hypoth- 
esis formulated in [30] that long term economic growth can be achieved by 
balancing key macroeconomic factors. In this example a five-year average of 
real GDP (gross domestic product) growth rate {Ggr) is used as a measure 
of long term economic growth, and the effect of the following macroeconomic 
indicators as independent variables is analyzed: inflation rate (/), balance of 
current account {BCA), and public debt {PD). An additional independent 
variable, GDP per capita (Gpp) , was used as a control variable to distinguish 
between different levels of economic development in the countries. 

Each point in our model thus represents a point in with coordinates the 
values of J, BCA and PD for one of the European countries in a given year 
from 1998 to 2004, and the average value of Gpp over this and the next four 
years. The average value of Ggr over this and the next four years is the depen- 
dent variable. Altogether 143 data points representing the countries included 
in different years were available. 

A Delaunay triangulation on the points in was constructed, producing 
altogether 17451 simplices of dimension up to 4. The discrete vector field was 
reconstructed from the data points using the algorithm of [12] producing 2 
critical cells in dimension 4 (maxima), 13 in dimension 3, 24 in dimension 2, 16 
in dimension 1 and 4 critical cells in dimension 0. On a laptop, the descending 
and ascending disks are reconstructed immediately, once the triangulation 
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and the discrete vector field are given. The reconstruction of the triangulation 
is also almost immediate, and the discrete gradient vector field is built in 
less than a minute. The graph, connecting the critical cells with respect to 
incidence in the boundaries of the corresponding descending regions is given 
in figure 14.41 



Figure 12. The qualitative graph connecting the critical cells in the 4-dimensional 
macroeconomic data set. 

For each country, data for the current year determines a point in the domain, 
and the model can be used to predict its average economic growth in the next 
five years. The V^-paths encode the changes in the parameter values which 
ensure constant growth of the GDP growth rate. Finally, the Morse-Smale 
decomposition can be used to construct a path, similar to the one in example 
14.31 which ends in a preferred position involving minimal variation in the GDP 
growth. 

The two maxima correspond to Ireland in 1999 (mi) and Latvia in 2003 (m2). 
According expert knowledge from economics, mi represents the preferred max- 
imum (since the high growth in Latvia is partly due to its post transition stage 
from a state planned economy, and involves high risks in the form of very high 
inflation, public debt and unemployment). So, for a country which is currently 
in the descending disc of 1712 (for example Slovenia), a path leading from its 
current position to a common 3-saddle, and from there along the gradient 
vector field to mi encodes the recommended strategy. 
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